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The objective of this effort is to develop an efficient and accurate computational 
methodology to predict both detailed and global thermo-fluid environments of a single flow 
element in a hypothetical solid-core nuclear thermal thrust chamber assembly. Several 
numerical and multi-physics thermo-fluid models, such as chemical reactions, turbulence, 
conjugate heat transfer, porosity, and power generation, were incorporated into an 
unstructured-grid, pressure-based computational fluid dynamics solver. The numerical 
simulations of a single flow element provide a detailed thermo-fluid environment for thermal 
stress estimation and insight for possible occurrence of mid-section corrosion. In addition, 
detailed conjugate heat transfer simulations were employed to develop the porosity models 
for efficient pressure drop and thermal load calculations. 

I. Introduction 

T he nuclear thermal rocket is one of the candidate propulsion systems for future space exploration including 
traveling to Mars and other planets of the solar system. Nuclear thermal propulsion can provide a much higher 
specific impulse than the best chemical propulsion available today. A basic nuclear propulsion system consists of 
one or several nuclear reactors that heat hydrogen propellant to high temperatures and then allow the heated 
hydrogen and its reacting product to flow through a nozzle to produce thrust. In the 1970's, a solid-core design * 1 * for 
the nuclear reactor was developed and tested under a nuclear rocket program called Rover/NERVA. The study 
showed that the solid-core reactor is a feasible concept to produce specific impulses exceeding 850 sec. The solid- 
core reactor resembles a heat exchanger, which consists of hundreds of heat generating solid flow elements, each 
flow element containing tens of flow channels through which the hydrogen propellant absorbs heat before entering a 
nozzle to generate thrust. To achieve maximum efficiency, the reactor often operates at very high temperature and 
power density. However, the results of Rover/NERVA tests indicated that if the solid fuel temperature exceeds 
certain values, the flow element may fail due to a phenomenon called mid-section corrosion 1 * 3 , which imposes real 
challenges to the integrity of the flow element. This limit imposes a constraint on the engine performance. The mid- 
section corrosion refers to a crack in the coating layer between the solid fuel and hydrogen flow, whichis designed to 
protect the solid fuel from chemical attack by the hot hydrogen. Mid-section corrosion can lead to an excessive mass 
loss of the flow element material in that region. The mid-section corrosion was suspected to be caused by different 
thermal expansions between the flow element and its coating, where the disparity can be large if strong thermal 
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gradients occur in that region, and where there is a change of solid thermal property due to irradiation. 4 Another 
speculation was that the flow was choked in the long flow channels of the flow element. 

It is extremely difficult and expensive to measure the detailed thermal-fluid 
environment within the entire flow element. The objective of this effort is 
therefore to develop an efficient and accurate multiphysics thermal-fluid 
computational methodology to predict environments for a single flow element, 
similar to those in the Small Engine. 1 The computational methodology was based 
on an existing Unstructured-grid Navier-Stokes Internal-external computational 
fluid dynamics Code (UNIC 5 ' 7 ). Conjugate heat transfer (CHT) formulations for 
coupling fluid dynamics and conductive heat transfer in solids and for flow and 
conductive heat transfer in porous media were developed and tested. The UNIC 
code has been well validated and employed to simulate a great variety of 
engineering problems ranging from internal to external flows, incompressible to 
compressible flows, single-phase to multi-phase flows, and inert to reacting flows. 

Two groups of detailed analyses were conducted for a NERVA-type 19-channel 
flow element, as shown in Figure 1. The first group simulated a full-length 19- 
channel flow element with two different power generation distributions to 
investigate the occurrence of the mid-section corrosion problem and potential 
flow chocking in the flow channel. These two distributions include a cosine 
function in the axial direction with a clipped cosine function in the radial 
direction, and a clipped cosine function in the axial direction with a clipped cosine 
function in the radial direction. The second group simulated a one-eighth-length 
19-channel flow element with different side-wall temperature and flow channel 
diameters to obtain data for calibrating the porosity model to be used in the global 
analysis of the entire thrust chamber. The porosity model is employed to 
represent the effect of the friction loss and heat transfer as the working fluid 
flowing through the flow element. The use of the porosity model enables an Figurel. Geometry of a 19- 

efficient simulation of the entire thrust chamber and evaluation of its performance channel flow element 
during the design cycle. 

In order to support both the detailed and global analyses, several physical submodels were either improved or added 
into the UNIC code 8 " 9 . These code improvements include an anisotropic drag and heat transfer porosity model to 
account for the directional effect of the flow channel, a source term in the energy equation to model the power 
generation of the solid core with user specified distributions, validation of the conjugate heat transfer capability with a 
well-known thermal analysis code, SINDA 10 . The porosity model development for the flow element and detailed 
conjugate heat transfer simulations of a powered flow element are reported herein. 

II. Numerical Methodology 

A. Computational Fluid Dynamics 

The employed CFD solver, UNIC, solves a set of Reynolds-averaged governing equations (continuity, Navier- 
Stokes, energy, species mass fraction, etc.) to satisfy the conservation laws for a turbulent flow of interest. The set of 
governing equation can be written in Cartesian tensor form: 
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where /> is the fluid density p is the pressure, Vj = (w, v, w) stands for the velocity components in >>-, and z- 
coordinates respectively, //, and h are the total and static enthalpies, k is the turbulence kinetic energy, P k and e are 
the production and dissipation rates of turbulence, a, and S, are the mass fraction and production/destruction rate of 
/- th species, Q r is the radiative heat flux, Sy and S h are the source/sink terms of the momentum and energy equations, 
p and p, are the fluid and eddy viscosity, t p represents the sum of the viscous and Reynolds stresses, Pr and Pr, are 
the Prandtl and turbulent Prandtl numbers. Sc and Sc, are Schmidt and turbulent Schmidt numbers, C h C 2 , C 3 , cr*, and 
c c are turbulence modeling constants. Detailed expressions for the k-e models and wall functions can be found in 
Ref. 11. An extended k-e turbulence model 12 was used to describe the turbulent flow. A modified wall function 
approach 13 14 was employed to provide wall boundary layer solutions that are less sensitive to the near-wall grid 
spacing. 

A predictor and multi-corrector pressure-based solution algorithm 15, 16 was employed in the UNIC code to 
couple the set of governing equations such that both compressible and incompressible flows can be solved in a 
unified framework without using ad-hoc artificial compressibility and/or a pre-conditioning method. The employed 
predictor-corrector solution method 5 is based on modified pressure-velocity coupling approach of the SIMPLE- 
type 16 algorithm which includes the compressibility effects and is applicable to flows at all speeds. In order to 
handle problems with complex geometries, the UNIC code employs a cell-centered unstructured finite volume 
method 6, 7 to solve for the governing equations in the curvilinear coordinates, in which the primary variables are the 
Cartesian velocity components, pressure, total enthalpy, turbulence kinetic energy, turbulence dissipation and mass 
fractions of chemical species. 

The inviscid flux is evaluated through the values at the upwind cell and a linear reconstruction procedure to 
achieve second order accuracy. A multi-dimensional linear reconstruction approach by Barth and Jespersen 1 was 
used in the cell reconstruction to achieve higher-order accuracy for the convection terms. A second-order central- 
difference scheme was employed to discretize the diffusion fluxes and source terms. A dual-time sub-iteration 
method is employed for time-accurate time-marching computations. A pressure damping term, Rhie and Chow 18 , is 
applied to the evaluation of mass flux at the cell interface to avoid the even-odd decoupling of velocity and pressure 
fields. All the discretized governing equations are solved using the preconditioned Bi-CGSTAB 19 matrix solver, 
except the pressure-correction equation which has an option to be solved using GMRES 20 matrix solver when the 
matrix is ill-conditioned. An algebraic multi-grid (AMG) solver 21 was included such that users can activate to 
improve the convergence if desired. In order to efficiently simulate problems involving large numbers of meshes, 
the UNIC code employed parallel computing with domain decomposition, where exchange of data between 
processors is done by using MPI 22 . Domain decomposition (partitioning the computational domain into several sub- 
domains handled by different computer processors) can be accomplished by using METIS 2 ' or a native partitioning 
routine in the UNIC code. 


B. Conjugate Heat Transfer (CHT): 

The framework of CHT is to solve the heat transfer in the fluid flow and the heat conduction in the solid in a 
coupled manner. The governing equation describing the heat transfer in the fluid flow is shown in Eq. (3), where the 
heat conduction in the solid can be written as: 
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where Q v and Q s represent source terms from volumetric and boundary contributions, respectively, /rand C p denote 
the thermal conductivity and capacity of the solid material, respectively. In the case of simulating solid core with 
power generation, Q v depends on power generation distributions function employed. In the present study, three 
power distribution functions were examined: (i) pure sine function in the axial direction with pure cosine function in 
the radial direction; (ii) pure sine function in the axial direction with clipped cosine function in the radial direction; 
(iii) Clipped sine function in the axial direction with clipped cosine function in the radial direction. The temperature 
value at the fluid-solid interface is obtained by enforcing the heat flux continuity condition, i.e. Q s = -q w where q w is 
the heat flux from the fluid to the solid calculated by solving Eq. (3) including the turbulence effect if applicable. In 
order to achieve numerical stability and enforce the heat flux continuity condition, an implicit treatment of the 
temperature at the fluid-solid interface is employed. In this approach, Eq. (8) can be discretized as 



where the superscripts rt + 1 and n denote the values at the next and current time levels, respectively. T h and T t are the 
temperatures at the fluid-solid interface and at the center of the solid cell next to the fluid-solid interface in respect. 
y n represents the normal distance from the interface to the center of the solid cell next to the fluid-solid interface. 
The Vi factor on the left hand side of the above equation is because only half of the solid cell is involved in the 
control volume. It can be seen that this scheme can be applied to both transient and steady state simulations. For 
steady-state simulations, an acceleration factor can be used to improve convergence of heat conduction in the solid. 
Implementation of the implicit treatment has been validated by comparing with the SINDA code. 

C. Porosity Model: 

In the present study, a porosity model was developed to represent the momentum and energy transport through 
an assembly of flow pipes and the heat conduction through the solid material within a flow element. Hence, the 
porosity model will include separate temperatures and thermal conductivities for both the solid material and fluid 
flow. The momentum and energy equations for the fluid flow are the same as Eqs. 2 & 3, except the source terms 
will be modified to account for the extra friction loss and heat transfer. The source term of the momentum equation 
(Eq. 2) can be expressed as 
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where D is the drag force modeled by the porosity model, E, is the porosity factor for the porous region, V, is the 
superficial flow velocity through the porous region, d is diameter of the flow channel, V is the total volume of the 
porous region, Cp is the drag coefficient, and f D is a modeling constant that will be tuned for the geometry of interest 
in the present study by comparing solutions of the porosity model and the detailed CHT model. An empirical 
correlation of the friction factor (Q) for the flow through a pipe, known as the Blasius formula 24 , is used and 
expressed as follows. 


C, =0.0791 Re d ° 25 where 
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The source term for the energy equation (Eq. 3) can be calculated as 

S =i- = A_^_c Af where AT = T,-T 0D 

* 2 Pr^ 3 gd / 4 p ' * 

where Q is the heat sink/source due to the fluid-solid interaction in the porous region, Pr is the Prandtl number of the 
fluid, f h is aa modeling constant will be tuned by comparing solutions of the porosity model and the detailed CHT 
model, and T s and T g are the temperatures of the solid and fluid at the same location, respectively. In addition, the 
source term of the heat conduction equation (Eq. 8) in the porous region can be calculated as 
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III. Numerical Meshes 

For the numerical simulations with the CHT model, two mesh systems need to be constructed. One of them is 
employed to model the entire computational domain including both the solid grain and the flow channel. The other 
one is used by the code to identify the region of solid material, and thus it can be very coarse. Hence, there is no 
need for a special grid generator to construct multi-block meshes for separating the solid and the fluid region, nor the 
identification for the fluid-solid interface. In the present study, various flow element geometries were simulated, 
such as different flow channel diameter, with and without coating layer in the solid core, and different flow channel 
length. Hence, a geometry/grid template, as shown in Figure 2, was developed based on MiniCAD 25 ' 2 to expedite 
the geometry and grid generation process by allowing the user to interactively change I ) the ratio of the gap between 
flow channels to the diameter of flow channel, 2) the ratio of the size of the flow element to the diameter of the flow 
channel, 3) thickness of the coating layer, and 4) the length of the flow element. The benefit of using MiniCAD is 
the ability to have graphical feedback while building the geometry and grid, and to easily move dimensional values 
from the geometry into the template to make it easy to create grids with modified values. In the present study, a 60° 
pi-section of a single flow element, which includes 3 and 1/6 flow channels, was simulated. 

MiniCAD is a graphical user interface (GUI) built on the Geometry Grid Toolkit (GGTK). The geometry that 
defines the dimensions of the solid grain was also exported as a discrete surface, which can be used by an in-house 
grid generator to construct the hybrid mesh for the entire flow element. The process of building grid in MiniCAD 
involves three basic steps. Step one is to create or import the geometry. Step two is to pass the geometry to a grid 
topology to ensure that the mesh built will be geometrically watertight. Step three is to build the mesh based on the 
grid topology. Based on our investigation, we found the most efficient structured grid topology to model the solid 
grain is an O-type grid. In this case, the solid grain was decomposed into multiple blocks (5 blocks for cases without 
a coating layer, and 1 0 blocks for those with a coating layer), where the O-grids of each block were used to model 
the solid grain and coating layer (if present) surrounding each flow channel as shown in Figure 2. A separate block 
is needed to identify the coating layer because it has different thermal properties from those of the solid grain. The 
boundary of the O-grid geometry for each flow channel was built as a planar trimmed surface. This geometry was 
converted to the topology structure, and then a 2-D grid was built on the faces. Since the grid for the solid grain is 
only used to identify the region of solid grain and is geometrically similar in the axial direction, it was only necessary 
to create the 3-D grid by extruding the 2-D grid in the axial direction, which can be seen in Figure 3. 



Figure 2. Layout of the developed geometry/grid template (flow channel diameter = 2.05 mm) 
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Figure 3. Structured grid system for the solid grain and coating layer of a flow element 

To accurately resolve the momentum and energy transport near the solid-fluid interface, a hybrid mesh topology 
was employed to generate the mesh for the entire flow element. The procedure used by our in-house grid generator 
is to construct the surface mesh before creating the 3-D volume mesh. For surface mesh generation, a direct 
advancing front method is used 0, °. A discrete model imported from MiniCAD is used as a background mesh, on 
which a new surface mesh suitable for computational simulations is created. Geometrical features are extracted 
based on a folding angle at each edge. A user specifies node distributions on the geometrical features, which form an 
initial front for the advancing front method. The surface triangulation is performed in the physical 3D space in order 
to check the quality of triangles easily. 

For volume mesh generation, two approaches were used. One is a method to create regular hybrid meshes based 
on an advancing layer method. In this method, a multiple marching direction approach is employed to create better 
control volumes around sharp convex comers 0 . Figure 4 shows a hybrid mesh inside a flow element, which has 
335,207 nodes, 688,027 tetrahedra, 21,476 prisms, 33,655 pyramids and 167,259 hexahedra. Quadrangles are used 
for the pipe surfaces. Multiple marching directions are defined at pipe inlets, which enable smooth transition of the 
mesh size there. 



Figure 4. Regular hybrid meshes: (a) surface mesh; (b) cross-section of hybrid volume mesh 


6 

American Institute of Aeronautics and Astronautics 


The other volume mesh generation method is an extrusion 
method based on a 2-D mesh to create volume meshes more 
easily. Every cross-section across the symmetry axis is the 
same. First, a 2-D mesh is created as shown in Figure 5. 

Second, it is extruded toward the symmetry axis, following the 
user-specified spacing. Third, boundary surfaces are 
numbered so that boundary conditions are set for CFD 
simulations. 

IV. Results and Discussion 

A. Calibration of porosity model: 

To tune the modeling constant employed in the porosity 
model, numerical experiments of simulating coolant flow 
through short-length flow elements with fixed boundary wall 
temperature (T*) using the CHT model were conducted to 
generate the anchoring data. The computational domain and 
boundary setup is shown in Figure 6. In this numerical experiment, the flow element has a width of 19.1 mm, length 
of 1 10 mm, and the inlet flow conditions are the same as those in Table 1, except the pressure is 34 atm. Even 
though the flow element was shortened from 890 mm to 1 10 mm to save the computational time, a length-diameter 
ratio (Lp/Dp) of 55 is enough to obtain a fully developed flow result for anchoring the porosity model. Two groups 
of numerical experiments were conducted. There are three cases in the first group, where the flow channel diameter 
is 2.05 mm and the wall temperatures are set to be 3000 K, 2000 K, and 1000 K, respectively. In the second group, 
the wall temperature was set to be 3000K, while three flow channel diameters (2.05 mm, 2.29 mm, and 2.54 mm) 
were modeled. A hybrid mesh system of 1.8 millions cells was used to simulate the coolant flow through a flow 
element using CHT. 



Figure 5. 2-D mesh for the extrusion method: 
solid part (yellow) and fluid part (green) 




(a) Cross-sectional view (b) Side view 

Figure 6. Sketch of the flow element geometry and its boundary conditions 



The numerical results of the first group (varying wall temperatures) and second group (varying flow channel 
diameter) follow the similarity rule, and thus only the flowfield of the case with a wall temperature of 3000 K and a 
flow channel diameter of 2.05 mm is plotted as shown in Figures 7-9. Figure 7 shows the pressure and temperature 
distributions at each boundary including the fluid-solid interface. The streamlines and velocity vectors of the fluid 
flow near the entrance and exit of the flow element are illustrated in Figure 8. Figure 9 shows the pressure 
distributions along the centerlines of the center and outer flow channels, and temperature distributions along 4 axial 
lines (Line #1: centerlines of the center flow channel; Line #3: the line between the center and outer flow channels; 
Line #5: centerlines of the outer flow channel; Line #7: the line between the outer flow channel and the edge of the 
flow element). As noted in Figure 9, a portion of temperatures along Lines #3 and #7 are the temperatures of the 
solid grain instead of the fluid flow. The temperature rises and pressure drops calculated from the numerical 
experiments will be used to tune the porosity model, and thus the results of other cases will be shown in the 
following comparison with the results of the porosity model. 
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TM 300 450 600 750 900 1050 1200 1350 1500 1650 1800 1950 2100 2250 2400 2550 2700 2850 3000 



P 33 84 33 86 33 88 33 9 33 92 33 94 33 96 33 98 34 34 02 34 04 34 06 34 08 34 1 

Figure 7. Temperature (top) and pressure (bottom) distributions for flow through a heated flow element 

with CHT (r H . = 3000 K, D p = 2.05 mm) 





Figure 8. Streamlines (left) and velocity vectors (right) near the entrance (top) and exit (bottom) of the flow 

element ( T w = 3000 K, D p = 2.05 mm) 
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Figure 9. Pressure (left) and temperature (right) distributions along the axial direction (T w = 3000 K, D p = 

2.05 mm) 

The cases simulated in the previous numerical experiments using CHT repeated using the porosity model. A 
hybrid mesh system with 23k cells was used to model the flow element, which is much coarser than that used in the 
simulation with CHT. The modeling constants of the employed porosity model, shown in Eqs 10 and 1 1, were tuned 
based on the cases of T w = 3000 K and D p = 2.05 mm, and were selected to be ^ = 18, an d f h = 0.22. The selection 



Axai distance along the centerfcne(m) Axal distance along the centerline (m) 
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of T w = 3000 K is because it is closest to the temperature calculated in the simulation of the full-length flow element 
with power generation and CHT. The calculated temperature and pressure distributions of this case are plotted as 
shown in Figure 10, which shows close resemblance to the result shown in Figure 7. With the same modeling 
constants, other cases were simulated with the porosity model. The results are compared to those obtained from the 
simulations using CHT, and are shown in Figures 11 and 12. Figure 11 indicates that the effect of the boundary 
temperature on pressure gradients calculated by the CHT model is more pronounced than that calculated by the 
porosity model. However, temperature distributions predicted by the porosity model agree very well with those by 
the CHT model. Variation of pressure gradients for different flow channel sizes calculated by the CHT and porosity 
models show minor discrepancy as exhibited in Figure 12. On the contrary, the result of the porosity' model displays 
a stronger influence of the flow channel size on temperature distributions than that by the CHT model. The cause for 
the discrepancy between the porosity and CHT models may be attributed to several factors. Since the local heat 
transfer and drag coefficients are highly dependent on the local flow Reynolds numbers, the value of the fluid 
viscosity need to be accurately accounted for. In the present study, the fluid viscosity of air, instead of hydrogen, 
was used, which can be one of the causes of the discrepancy. The accuracy of the correlation between the fluid 
viscosity and temperature employed in the code also needs to be further examined. The Blasius formula for the 
friction loss coefficient, developed based on the isothermal flow, may cause some errors of the non-isothermal fluid 
flow in the present study, and thus requires further refinement. The result shows the need for further investigation of 
the employed empirical correlations of the porosity model such that it can be applied to a wider variation of flow 
conditions and flow element geometries. Despite this discrepancy, if the operating conditions and geometry are 
close to the tuning conditions (maximum T w = 3000 K, D p = 2.05 mm), it is appropriate to use the tuned porosity 
model to predict the global performance of nuclear thermal thrust chamber assembly. 



TM 300 450 600 750 900 1050 1200 1350 1500 1650 1800 1950 2100 2250 2400 2550 2700 2850 3000 



P 33 89 33 91 33 93 33 95 33 97 33 99 34 01 34 03 34 05 34.07 34 09 34 11 34 13 34 15 

Figure 10. Temperature (top) and pressure (bottom) distributions for flow through a heated flow element 

with porosity model (T w = 3000 K, D r = 2.05 mm) 




temperatures (D p = 2.05 mm) 
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Figure 12. Comparisons of pressure gradient (left) and temperature (right) distributions with various flow 

channel diameters (T w = 3000 K) 

B. Full-length single flow element with power generation: 

To investigate the mid-section corrosion problem, a coolant flow (gaseous hydrogen) though the full-length, 
innermost flow element (the largest energy output) with a given power generation source was simulated to 
understand the thermal and fluid environment with the flow element. The layout of the single flow element 
simulated and its geometric definitions are similar to those illustrated in Figure 6, where only 1/6 of the flow element 
(shaded area in Figure 6) is considered in the numerical simulation due to geometrical symmetry. The only 
difference is that instead of an isothermal condition (at T w \ an adiabatic boundary condition was imposed at the edge 
of the flow element. A hybrid mesh system, as shown in Figures 4-5, with 4.5 millions cells was used to model the 
flow element, the upstream and downstream sections. The values of the geometric dimension of the flow element, 
and the inflow conditions were obtained from Ref. 31, and are listed in Table 1. A coating layer (ZrC composite) 
with a thickness of 0.127 mm was added to each flow channels. In the numerical simulation, the thermal properties 
of the flow element (UC-C-ZrC composite) and the coating layer (ZrC composite) vary with temperature and can be 
found in Ref. 2, where the values of those properties at a temperature of 300 K are listed in Table 2 as an example. 
The thermal properties (UC-C-ZrC composite 2 ) are listed in Table 2. It is noted that in the actual simulation the 
thermal properties of UC-C-ZrC and ZrC composites vary with temperature, where the values listed in Table 2 are 
the properties at a temperature of 300 K. An entrance section of 1 element width (W p ) was added to the upstream of 
the flow element, and an exit section of 2 W p was added to the downstream of the flow element. For the power 
generation from the solid fuel grain, a clipped cosine distribution in the radial direction (maximum power at the 
center of the innermost flow element, and zero power at the edge of the outermost flow element) was specified. Two 
distribution curves of power generation ($) in the axial direction were investigated. Case #1 has a pure sine 
distribution and Case #2 has a clipped sine distribution, and both are shown in the following 

Case #1 : </> = ^ sin(;rx / L p ) x :0 L p 

Case #2 : (f> = <t> mMX min[ 1, 1 .5 sin(^x / L p )] x : 0 -> L p 

where (f^ is the maximum power and L p is the length of the flow element. The total power generated is set to be 
491 kW, which is about 80% of the designed value 31 . The reason for selecting this value is because this power is 
high enough to cause flow chocking in the flow' element. Further increase of the power will lead to unrealistic 
numerical solutions because of the fixed mass inlet boundary condition and mass conservation exit boundary 
condition. This set of boundary conditions was used because only a 1/6 segment of one flow element was simulated 
and thus the nozzle (which fits multiple flow elements) can not included. Since each flow element has different 
energy input, the centermost flow element, subjected to highest heat addition, is simulated in the present study as the 
worst scenario. Furthermore, the purpose of this study is to investigate possible occurrence of flow chocking within 
the flow channel, what is the maximum flow Mach number within the flow channel at the designed power level is a 
secondary issue and will be a subject for the future study. The heat transfer between the fluid flow and the solid fuel 
is simulated using the CHT model. A single elementary reaction was employed to model the hydrogen dissociation 
process due to the high temperature. 
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Flow Element Geometry 

Inflow conditions (H 2 ) 

Element width 
(W P ) 

Channel 
diameter ( D p ) 

Element length 

(L P ) 

Pressure 

Temperature 

Mass flow rate 

19.1 mm 

2.05 mm 

890 mm 

39.05 atm 

279 K 

0.0147 kg/s/element 


Table 1. Geometry of a flow element and coolant flow inlet conditions 



Thermal conductivity 

Density 

Specific heat 

UC-C-ZrC composite 

91 W/m-K 

6142 kg/m 3 

468.4 W-s/kg-K 

ZrC composite 

19 W/m-K 

6531 kg/m 3 

368 W s/kg-K 


Table 2. Thermal properties of flow element and coating layer at 300 K 


The numerical results of power generation cases #1 and #2 are plotted as shown in Figures 13-18. Figure 13 
shows the pressure, temperature and Mach number contours at the symmetry boundaries and solid-fluid interfaces of 
the power generation case #1, while Figure 16 shows those features of the power generation case #2, respectively. 
As can be seen from the Mach number contours, the coolant flow is choked near the end of the flow channel. It is 
noted that these two plots are re-scaled due to extremely large aspect (length-to-width) ratios such that the contours 
are viewable and distinguishable. Figure 14 shows the pressure distributions along the centerlines of the center and 
outer flow channels, and temperature distributions along 7 axial lines of the power generation case #1. The location 
of each line is depicted in Figure 5 (Line #1: centerlines of the center flow channel; Line #2: the line along the mid- 
point of the coating layer around the center flow channel; Line #3: the line along the mid-point between the center 
and outer flow channels; Line #4: the line along the inner mid-point of the coating layer around the outer flow 
channel; Line #5: centerlines of the outer flow channel; Line #6: the line along the outer mid-point of the coating 
layer around the outer flow channel; Line #7: the line between the outer flow channel and the edge of the flow 
element). In the temperature distribution plot, the location of the maximum temperature gradient along each line is 
also marked with a symbol which has the same color as each line. It can be seen that Lines #2 and #4 have the same 
maximum temperature gradient locations, which are very close to that of Line #3. The axial temperature distribution 
is qualitatively similar to that of 1-D calculation. Some researchers think that the location of maximum temperature 
gradient is critically linked to the potential corrosion of the flow element caused by different thermal expansion 
coefficients between the coating layer and the fuel. Hence, temperature distributions in the radial direction (across 
the solid-fluid interface) at two maximum temperature gradient locations (of the coating layer and center flow 
channel) are also plotted as shown in Figure 15. These two axial locations can be seen in Figure 14. The radial 
temperature profiles reveal that substantial thermal gradients occur at the fluid-coating interface and the coating-fuel 
interface. This may be attributed to large difference of the thermal conductivities between the coating material and 
the fuel, as shown in Table 2. Though 1-D calculations are efficient and can provide estimation of axial temperature 
distributions and maximum thermal gradient locations, whether the mid-section corrosion will occur can only be 
speculated because no information about the temperature gradients across different materials and fluid is available 
for thermal stress and fracture analyses. This shows the need for conducting 3-D CFD simulations in order to obtain 
a detailed thermal-fluid environment in the flow element. It is noted that this is a steady-state result, and thus the 
thermal gradient could be much stronger during transient startup as pointed out by Wang et al 9 Figure 17 shows 
axial temperature distributions of the power generation case #2, which are similar to those of Figure 14 except that 
the locations of maximum thermal gradient are further upstream than those of case #1. This is expected because 
power generation case #2 has a clipped sine distribution, where the power rises up to its peak value faster in the axial 
direction. Temperature distributions in the radial direction at the same two maximum thermal gradient locations are 
plotted in Figure 18. Substantial thermal gradients at the fluid-coating interface and the coating-fuel interface also 
can be observed. However, the level of thermal gradient at the coating-fiiel interface is relatively smaller than that of 
the power generation case #1. This can be attributed to the characteristics of the clipped cosine distribution, which 
has a larger area subjected to the peak power source but has a smaller peak power (such that the overall energy 
inputs for both cases are the same). The exit temperatures of both cases are almost identical, and have a value about 
2660 K. 
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Radial distance (mm) Pressure (atm) 



P 7 10 13 16 19 22 25 28 31 34 37 40 43 46 



TM 300 450 600 750 900 1050 1200 1350 1500 1650 1800 1950 2100 2250 2400 2550 2700 2850 



M 01 02 03 04 05 06 0 7 08 09 1 1 1 

Figure 13. Pressure (top, in atm), temperature (middle, in K), and Mach number (bottom) contours at the 
boundary and solid-fluid interfaces of power generation case #1 
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igure 14. 
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Pressure (left) and temperature (right) distributions along the axial direction (power generation 

case #1) 






Figure 15. Radial temperature profiles at the maximum temperature gradient locations of coating layer 
(left) and center flow channel (right) for power generation case #1 
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Radial distance (mm) T) Pressure (atm) 



P 7 10 13 16 19 22 25 28 31 34 37 40 43 46 



TM 300 450 600 750 900 1050 1200 1350 1500 1650 1800 1950 2100 2250 2400 2550 2700 2850 



M 0 1 0 2 0 3 0 4 0 5 0.6 0.7 0 8 0 9 1 1.1 

Figure 16. Pressure (top, in atm), temperature (middle, in K), and Mach number (bottom) contours at the 
boundary and solid-fluid interfaces of power generation case #2 



igure 17. Pressure (left) and temperature (right) distributions along the axial direction (power generation 

case #2) 




Figure 18. Radial temperature profiles at the maximum temperature gradient locations of coating layer 
(left) and center flow channel (right) for power generation case #2 
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V. Conclusions 

In the present study, the developed CFD model (including conjugate heat transfer, power generation, and 
porosity model) has been successfully applied to simulate the thermal-fluid environment of a single flow element in a 
nuclear thermal thrust chamber. The use of the porosity model shows a promise of being an efficient numerical 
approach for simulating an assembly of a huge number of flow elements. However, the employed porosity model 
fails to match the numerical data of the CHT model at a wide range of operating and geometric conditions. Further 
study is required to improve the porosity model such that it can be employed to accurately simulate the global 
thermal-fluid environment of the flow element. A detailed numerical analysis of a powered flow element using the 
CHT model provides the insight of the thermal-fluid environment The numerical result indicates large thermal 
gradients at the fluid-coating and coating-fuel interfaces for the selected fuel and coating materials, which can 
potentially cause the mid-section corrosion problem. Moreover, the radial temperature profiles show that 3-D CFD 
simulations provide critical information which can not be obtained from 1-D calculations. Selection of different 
material may be able to reduce the level of thermal gradients. The results also reveal that the power distribution 
function can affect the location and level of maximum thermal gradient. Nevertheless, the numerical results of both 
power distributions indicate that even at 80% of the designed power level some flow elements possibly have flow 
chocking in the flow channel and some may not depending on the local heat source magnitude. The occurrence of 
flow chocking in some flow channels can cause non-uniform flow distributions among flow elements, and produce 
undesirable side load. Further study using 3-D CFD simulations of the full configuration of the flow element 
assembly with a downstream nozzle is needed in order to evaluate 1) the magnitudes of the maximum flow Mach 
number, 2) distributions of mass flow rates pressure drops among the flow elements, and 3) the magnitude of side 
load and its impact. Different flow element configurations should also be examined in the future to eliminate the 
flow chocking problem. 
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BACKGROUND 


• Nuclear thermal rocket is one of the candidate propulsion 
system for future space exploration 

• Study from Rover/NERVA program showed the solid-core 
reactor is a feasible concept to produce l sp over 850 sec 

• The solid-core reactor consists of hundreds of heat 
generating solid flow elements 

• Each flow element contains tens of flow channels through 
which H 2 is heated up before entering a nozzle 

• To achieve maximum efficiency: 

■ The reactor operates at very high temperature and power density 

■ Number and area of flow channels per flow element need to be 
optimized 


BACKG ROUND (Cont.) 


• A coating layer was designed to protect the solid fuel from 
chemical attack by hot H 2 

• Results of Rover/NERVA tests indicate: 

■Mid -section corrosion (crack in coating layer) occurs if the solid 
fuel temperature exceeds certain value, which leads to mass loss 
of the flow element 

■Mid -section corrosion was suspected to be caused by different 
thermal expansions between the flow element and coating layer 

■Dispa rity of thermal expansion can be large if large thermal 
gradients occur near the interface between the flow channel, 
coating layer and flow element 

■Possible flow choking in some of the flow channels, which will lead 
to non-uniform propellant flow distributions among flow elements 
and undesirable side load 



MOTIVATIONS 


• Performance of the solid-core reactor was previously 
estimated based on 1-D analysis, which lacks of detailed 
thermal-fluid environments to address the afore- 
mentioned problems 

• 3-D CFD simulations can provide high-fidelity thermal- 
fluid environments in the reactor, but need to 

■Imp rove computational efficiency 

■Couple calculation of heat transfer in the fluid flow with calculation 
of heat conduction in the solid fuel 

■Includ e realistic power generation model 

■Inco rporate correct thermal properties of the flow element and 
coating layer 

■Accou nt for dissociation of H 2 at high temperature 

■Accou nt for neutronic reactions 


OBJECTIVES 


• Develop an efficient, accurate CFD methodology to 
simulate thermal-fluid environments of a single flow 
element 

■ Improve the conjugate heat transfer model and implement the 
power generation model to accurately predict detailed thermal-fluid 
environment within flow element 

- Temperature, pressure, velocity distributions and possible choking of 
coolant flow 

- Thermal gradient in the solid and at the solid-fluid interface for possible 
corrosion/cracking 

■ Provide temperature distributions in the flow element for the 
thermal stress analysis 

■ Improve and tune the porosity to efficiently calculate global thermal- 
fluid phenomena: overall pressure loss & temperature rise of 
coolant flow 


Numerical Approach & Models 

• Unstructured-grid pressure-based method for solving Navier- 
Stokes or RANS equations 

• Cell center, finite volume approach 

• Compressible and incompressible flows with all speed regimes 

• Barth-Jespersen’s linear reconstruction approach 

• Conjugate Gradient, or GMRES Matrix Solvers 

• PISO scheme for unsteady flow simulations 

• Algebraic multi-grid (AMG) method to enhance convergence 

• Flexible mesh adaptation refinement 

• Parallel computing with MPI or PVM 

• Dynamic memory allocations 

• Domain decomposition with METIS 



Numerical Approach & Models (cont.) 

• Standard and extended k-s model for modeling turbulence 
effects 

• Finite-Rate and Equilibrium Chemistry Models 

• Conjugate heat transfer model for coupling heat convection 
and conduction 

• Heterogeneous spray model with Eulerian/Lagrangian 
framework and volume of fluid approach 

• Homogeneous spray model with Eulerian/Eulerian framework 
and real-fluid property model 

• Finite volume method for solving radiation transport equation 
with either gray gas or narrow band model 

• Porosity model for flow through porous media 

• Account for thermal non-equilibrium effect by solving 
vibrational energy equation 
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Geometry of Flow Element Assembly 
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Tuning of Porosity Model 



•T w : 1000, 2000, and 3000 K 
•D p \ 2.05, 2.29, and 2.54 mm 
•L p \ 110 mm, W p \ 19.1 mm 

• Flow element: (UC-C-ZrC composite) 

k= 22.5 W/m K, p= 13500 kg/m 3 , C p = 146.5 W s/kg-K 

• Inlet conditions: 

Pressure: 34 atm, Temperature: 279 K, Flow rate: 0.0147 kg/s 



Tuning of Porosit y Model 

Simulation with conjugate heat transfer model 


No. of cells: 1.7M, no. of grid points: 555K 






Tuning of Porosity Model 
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Tuning of Porosity Model 




Simulation with porosity model 


No. of cells: 22K, no. of grid points: 8.6K 
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Effect of Pore Size 










• No. of cells: 4.5M, no. of grid points: 98K 

• Geometry and inlet B.C. 


Flow Element Geometry 

Inflow conditions (H 2 ) 

Element 
width (W p ) 

Channel 
diameter (D p ) 

Element 
length (£ 0 ) 

Pressure 

Temperature 

Mass flow 
rate 

19.1 mm 

2.05 mm 

890 mm 

39.05 atm 

279 K 

0.0147 kg/s 


• Variable thermal properties for flow element and coating 
layer (listed are reference values at 300 K) 
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Analysis of Full-Si ze Flow Element 

• Case #1 : pure cosine profile for power generation 

Q = cos \n{x!L p - 0.5)] x : 0 > L /t 

• Case #2: clipped cosine profile for power generation 

^ = ^max min {U 1.5cos[;r(x/ L p -0.5)]] x:0-> L p 

•Total power generated: 491 kW (80% designed value) 
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Analysis of Full-Size Flow Element 
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Case #1 














Analysis of Full-Siz e Flow Element 

Case #1 




At max. temperature gradient of coating At max. temperature gradient of center 

layer flow channel 
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Analysis of Full-Size Flow Element 


Case #2 







Analysis of Full-Size Flow Element 
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Case #2 








Analysis o f Full-Size Flow Element 


Case #2 



At max. temperature gradient of coating At max. temperature gradient of center 

layer flow channel 





CONCLUSIONS 


• The employed CFD model has been successfully applied 
to simulate thermal-fluid environment of a single flow 
element 

• The use of porosity model can be an efficient numerical 
approach for simulating an assembly of a huge number of 
flow elements, but needs further improvement 

• The numerical result indicates large thermal gradients at 
the fluid-coating and coating-fuel interfaces for the 
selected fuel and coating materials 

• Power distribution function can affect the location and 
level of max. thermal gradient 

• Even at 80% of the designed power level some flow 
elements possibly have flow chocking in some flow 
channels 


RECOMMENDATIONS 

• A pie-section of the reactor with a downstream nozzle 
should be simulated to: 

■ Account for the multiple flow element effect 

■ Evaluate the magnitude of max. flow Mach no. in the flow channel 
with the full power generated 

■ Calculate distribution of the mass flow rate, temperature and 
pressure among flow elements 

• Conduct thermal stress analysis to examine the 
possibility of crack of the coating layer 

• Different configurations of the flow element (number and 
area) should be examined to optimize its efficiency 

• Neutronic reactions should be included to accurately 
model power generation 


